% Data analysis
clc
clear
load duration1quarters_TP06QRJfinR5G1000kMR7.mat;

% evolution of rating universe over time-----------------------------------

numIssuers = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, 1);

for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:, :) = durationCount(i,j,:,:);
        numIssuers(4*(i-1)+j,:) = sum(sum(x));
    end
end
numIssuersTS = timeseries(numIssuers(:,:), 1: length(numIssuers), 'name', 'Issuers Quarter');

figure
plot (numIssuersTS,'-*');
set(gca,'XTick',1:4:length(numIssuers));
labels = quaterlabels(1998, length(numIssuers));
labels = labels(1:4:length(numIssuers));
set(gca,'XTickLabel',labels);
%grid minor;
xlabel('Quarter');

% Average Rating Distribution----------------------------------------------
% Evolution of Rating Over Time

rateDist = zeros (RATE_LIST_LENGTH-1, 1);
rateDist06Before = zeros (RATE_LIST_LENGTH-1, 1);
rateOverTime = zeros( NUMBER_OF_YEAR * Quarter_LIST_LENGTH , RATE_LIST_LENGTH-1 );
for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:, :) = durationCount(i,j,:,:);
        rateDist = sum(x,2)+ rateDist;  % exclude withdrawn, how many obligors quarters totally.  
        if i==8 
            rateDist06Before = rateDist; % exclude withdrawn, how many obligors quarters totally before 2006q1.
        end          
        rateOverTime (4*(i-1)+j,:)= sum(x,2)';
    end
end

rateDist06After = rateDist - rateDist06Before; % exclude withdrawn, how many obligors quarters totally after 2006q1.

rateDistAvg = rateDist./(NUMBER_OF_YEAR*Quarter_LIST_LENGTH);
rateDistAvgBefore2006 = rateDist06Before./(8*Quarter_LIST_LENGTH);
rateDistAvgAfter2006 = rateDist06After./(4*Quarter_LIST_LENGTH);


figure;
bar (rateDistAvg);
set(gca,'XTickLabel',{'1.5','2','2.5','3','3.5','4','4.5' });
xlabel('Ratings');
ylabel('Obligors Quarter');
title('Average Rating Distribution');

figure;
bar (rateDistAvgBefore2006);
set(gca,'XTickLabel',{'1.5','2','2.5','3','3.5','4','4.5' });
xlabel('Ratings');
ylabel('Obligors Quarter');
title('Average Rating Distribution Before 2006Q1');

figure;
bar (rateDistAvgAfter2006);
set(gca,'XTickLabel',{'1.5','2','2.5','3','3.5','4','4.5' });
xlabel('Ratings');
ylabel('Obligors Quarter');
title('Average Rating Distribution After 2006Q1');

figure;
plot (rateOverTime);
set(gca,'XTick',1:4:length(numIssuers));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Obligors Quarter');
legend('1.5','2','2.5','3','3.5','4','4.5','Location','NorthWest');


% Default Freqency Evolution-----------------------------------------------
% Evolution of Default Frequency Over Time and Ratings

defaultFreq = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, 1);
defaultFreqOverRate = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, RATE_LIST_LENGTH-1);

for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        x(:, :) = durationCount(i,j,:,:);
        defaultFreq(4*(i-1)+j,:) = sum(x(:,RATE_LIST_LENGTH)) / sum(sum(x));
        defaultFreqOverRate(4*(i-1)+j,:) = ( x(:, RATE_LIST_LENGTH) ./ sum(x, 2) )';
    end
end

defaultFreqTS = timeseries(defaultFreq(:,:), 1: length(defaultFreq), 'name', 'Default Freqency');

figure;
plot (defaultFreqTS,'-*');
set(gca,'XTick',1:4:length(defaultFreq));
set(gca,'XTickLabel',labels);
grid minor;
xlabel('Quarter');

figure;
plot (defaultFreqOverRate);
set(gca,'XTick',1:4:length(defaultFreq));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Default Freqency');
legend('1.5','2','2.5','3','3.5','4','4.5','Location','NorthWest');

% Unconditional Transition Matrix------------------------------------------


stdProbQ = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
stdProbAvgQ = zeros (RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
mobilityNorm = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
mobilityNormRate = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
SVD = zeros(NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
SVDRate = zeros(NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
xx = zeros (RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);


for i = 1: NUMBER_OF_YEAR
    for j = 1: Quarter_LIST_LENGTH
        
        x(:, :) = durationCount(i,j,:,:);
        xx = xx + x;  
        
        for d = 1: RATE_LIST_LENGTH-1            
            stdProbQ (4*(i-1)+j, d, :) = x(d,:)./ sum(x(d,:),2);            
        end
                
    end
end

for d = 1: RATE_LIST_LENGTH-1
        stdProbAvgQ (d, :) = xx(d,:)./ sum(xx(d,:),2);        
end


% Metrics Calculation
% SVD (define singular value of matrix as mobility norm, here we use average of singular
% value instead of largest one which is argued in Jafry and Schuerman)
% D3 (risk-adjusted indices)

% NOTICE !! Here all the metrics is distance to average matrix (all sample)
L1 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
L2 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
WAD = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
NAD = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D3 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D1 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D2 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D4 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D1sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D2sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D3sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
D4sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);


rateMatrix = stdProbQ (:, :, 1:RATE_LIST_LENGTH-1);
rateMatrixAvg = stdProbAvgQ (:, 1:RATE_LIST_LENGTH-1);
a = mean(svd(stdProbAvgQ));
aRate = mean(svd(rateMatrixAvg));

for i = 1: size(stdProbQ, 1)
    aa(:,:) = stdProbQ (i, :, :);
    aaRate (:,:) = rateMatrix (i, :, :);
    mobilityNorm (i,1) = mean(svd(aa));
    mobilityNormRate(i,1) = mean(svd(aaRate));
    SVD(i,1) = mobilityNorm (i,1) - a; %SVD
    SVDRate(i,1) = mobilityNormRate (i,1) - aRate;
    
    for j = 1 : RATE_LIST_LENGTH-1
        for k = 1 : RATE_LIST_LENGTH
            
            L1(i) = L1(i) + abs(aa(j,k) - stdProbAvgQ(j,k)); % L1
            L2(i) = L2(i) + (aa(j, k) - stdProbAvgQ(j, k))^2;
            WAD(i)= WAD(i) +  abs(aa(j,k) - stdProbAvgQ(j,k))* aa(j,k); %WAD
            
            
            if aa(j,k)~=0
                NAD(i) =  NAD (i) + abs(aa(j,k) - stdProbAvgQ(j,k))/aa(j,k); %NAD
                mD2 = (j - k)* (aa(j, k) - stdProbAvgQ(j, k)) / aa(j, k);
                mD4 = (j - k)* sign(aa(j, k) - stdProbAvgQ(j, k))*(aa(j, k) - stdProbAvgQ(j, k))^2 / aa(j, k);
            end
            
            
            
            mD3 = (j - k)* sign(aa(j, k) - stdProbAvgQ(j, k))*(aa(j, k) - stdProbAvgQ(j, k))^2;
            mD1 = (j - k)* (aa(j, k) - stdProbAvgQ(j, k));
            
            if k~= RATE_LIST_LENGTH %D3
                D3(i) = D3(i) + mD3;
                D1(i) = D1(i) + mD1;
                D2(i) = D2(i) + mD2;
                D4(i) = D4(i) + mD4;
                D1sqr(i)= D1sqr(i)+ mD1;
                D2sqr(i)= D1sqr(i)+ mD2;
                D3sqr(i)= D1sqr(i)+ mD3;
                D4sqr(i)= D1sqr(i)+ mD4;
                
            else
                D3(i) = D3(i) + mD3*RATE_LIST_LENGTH;
                D1(i) = D1(i) + mD1*RATE_LIST_LENGTH;
                D2(i) = D2(i) + mD2*RATE_LIST_LENGTH;
                D4(i) = D4(i) + mD4*RATE_LIST_LENGTH;                
                D1sqr(i) = D1sqr(i)+ mD1*RATE_LIST_LENGTH^2;
                D2sqr(i) = D2sqr(i)+ mD2*RATE_LIST_LENGTH^2;
                D3sqr(i) = D3sqr(i)+ mD3*RATE_LIST_LENGTH^2;
                D4sqr(i) = D4sqr(i)+ mD4*RATE_LIST_LENGTH^2;
            end
            
            
        end
    end
    L2(i) = sqrt(L2(i)); %L2
end


% 
% figure;
% plot (mobilityNorm);
% set(gca,'XTick',1:4:length(mobilityNorm));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Mobility Norm SVD');

figure;
subplot(4,2,1); plot (SVD);
set(gca,'XTick',1:4:length(mobilityNorm));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('SVD');

subplot(4,2,2);plot (D3);
set(gca,'XTick',1:4:length(D3));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Risk-adjusted Metric: D');

subplot(4,2,3);plot (L1);
set(gca,'XTick',1:4:length(L1));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('L1');

subplot(4,2,4);plot (L2);
set(gca,'XTick',1:4:length(L2));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('L2');

subplot(4,2,5);plot (WAD);
set(gca,'XTick',1:4:length(WAD));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('WAD');

subplot(4,2,6);plot (NAD);
set(gca,'XTick',1:4:length(NAD));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('NAD');

subplot(4,2,7);plot (D1);
set(gca,'XTick',1:4:length(WAD));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('D1');

subplot(4,2,8);plot (D1sqr);
set(gca,'XTick',1:4:length(NAD));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('D1sqr');



% Diagnoal

diagProbQ = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, RATE_LIST_LENGTH-1);

for q = 1 : NUMBER_OF_YEAR * Quarter_LIST_LENGTH
    for d = 1: RATE_LIST_LENGTH-1        
        diagProbQ (q,d) = stdProbQ (q, d, d);        
    end
end

figure;
plot (diagProbQ);
set(gca,'XTick',1:4:length(diagProbQ));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Non-Transition Frequency');
legend('1.5','2','2.5','3','3.5','4','4.5','Location','NorthWest');


% Distribution of PD and inverse normal CDF of PD--------------------------
% Distribution of SVD and stdSVD
% Distribution of D3
% Distribution of L1
% Distribution of WAD, NAD
% 

% Compare the actual distribution with normal distribution.

defFreqInverse = norminv(defaultFreq);
stdSVD = -zscore(SVD);


figure
subplot (1, 2, 1);
probplot(defaultFreq);
legend('Normal','PD','Location','NW')
subplot (1, 2, 2);
probplot(defFreqInverse);
legend('Normal','Inverse Normal CDF of PD','Location','NW')

figure

subplot (4, 2, 1);
probplot(L1);
legend('Normal','L1','Location','NW');
subplot (4, 2, 2);
probplot(L2);
legend('Normal','L2','Location','NW')
subplot (4, 2, 3);
probplot(WAD);
legend('Normal','WAD','Location','NW')
subplot (4, 2, 4);
probplot(NAD);
legend('Normal','NAD','Location','NW')
subplot (4, 2, 5);
probplot(SVD);
legend('Normal','SVD','Location','NW');
subplot (4, 2, 6);
probplot(D3);
legend('Normal','D3','Location','NW')
subplot (4, 2, 7);
probplot(D1);
legend('Normal','D1','Location','NW');
subplot (4, 2, 8);
probplot(D1sqr);
legend('Normal','D1sqr','Location','NW')

% probplot(defaultFreq)
% p = mle(defaultFreq,'dist','tlo');
% t = @(defaultFreq,mu,sig,df)cdf('tlocationscale',defaultFreq,mu,sig,df);
% h = probplot(gca,t,p);
% set(h,'color','r','linestyle','-')


% chi-square goodness-of-fit test for testing normal distribution of PD and
% inverse normal CDF of PD
[h(1),p(1)] = chi2gof(defaultFreq);
[h(2),p(2)] = chi2gof(defFreqInverse);
[h(3),p(3)] = chi2gof(L1);
[h(4),p(4)] = chi2gof(L2);
[h(5),p(5)] = chi2gof(WAD);
[h(6),p(6)] = chi2gof(NAD);
[h(7),p(7)] = chi2gof(SVD);
[h(8),p(8)] = chi2gof(D3);
[h(9),p(9)] = chi2gof(D1);
[h(10),p(10)] = chi2gof(D1sqr);


% draw graph for SVD, D1-D4 , D1_n^2 - D4_n^2, (L1,WAD,NAD) distance to average TM
t = 1:1:NUMBER_OF_YEAR * Quarter_LIST_LENGTH;

figure;
plot(t,0); hold on % add recession figure on previous graph
ha = area([41 47], [0.03 0.03]);
set(ha,'BaseValue',-0.05);
set(ha,'FaceColor',[0.83 0.82 0.78]);
set(ha,'LineStyle','none')
grid on
set(gca,'Layer','top'); hold on % add metrics on previous graph
hb = plot (SVD,'-+');
set(gca,'XTick',1:4:length(SVD));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Mobility Norm SVD');
grid on

figure;
plot(t,0); hold on % add recession figure on previous graph
ha = area([41 47], [5 5]);
set(ha,'BaseValue',-4);
set(ha,'FaceColor',[0.83 0.82 0.78]);
set(ha,'LineStyle','none')
grid on
set(gca,'Layer','top'); hold on % add metrics on previous graph
hb = plot(t,zscore(D1),'-+',t,zscore(D2),'-o',t,zscore(D3),'-*',t,zscore(D4),'-x');
legend(hb,'D1','D2','D3','D4','Location','NorthWest');
set(gca,'XTick',1:4:length(numIssuers));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Risk-adjusted Difference Metrics D1-D4');
grid on


figure;
plot(t,0); hold on % add figure on previous graph
ha = area([41 47], [4 4]);
set(ha,'BaseValue',-3);
set(ha,'FaceColor',[0.83 0.82 0.78]);
set(ha,'LineStyle','none')
grid on
set(gca,'Layer','top'); hold on
hb = plot(t,zscore(D1sqr),'-+',t,zscore(D2sqr),'-o',t,zscore(D3sqr),'-*',t,zscore(D4sqr),'-x');
legend(hb,'D1(n^2)','D2(n^2)','D3(n^2)','D4(n^2)','Location','NorthWest');
set(gca,'XTick',1:4:length(numIssuers));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Risk-adjusted Difference Metrics D1_n^2-D4_n^2');
grid on

figure;
plot(t,0); hold on % add figure on previous graph
ha = area([41 47], [4 4]);
set(ha,'BaseValue',-2);
set(ha,'FaceColor',[0.83 0.82 0.78]);
set(ha,'LineStyle','none')
grid on
set(gca,'Layer','top'); hold on
hb = plot(t,zscore(L1),'-+',t,zscore(WAD),'-o',t,zscore(NAD),'-*');
legend(hb,'L1','WAD','NAD','Location','NorthWest');
set(gca,'XTick',1:4:length(numIssuers));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
ylabel('Risk-adjusted Difference Metrics L1,WAD,NAD');
grid on


% draw the graph of cell prob along time
% We did the same cell as paper Xing (2010)
% ------------------
% 1 AAA 1.5
% 2 AA  2
% 3 A   2.5
% 4 BBB 3
% 5 BB  3.5
% 6 B   4
% 7 C   4.5
% 8 D   D
%-------------------
% the cell to check include: 1-1, 5-5, 1-8, 5-8, 2-6, 3-4

figure;
subplot(3,2,5); plot (stdProbQ(:,2,2));
legend('Rating 2 to 2 ');
title('Observed Transition Prob along Time',... 
  'FontWeight','bold')
set(gca,'XTick',1:4:length(numIssuers));
labels = quaterlabels(1998, length(numIssuers));
labels = labels(1:4:length(numIssuers));
set(gca,'XTickLabel',labels);
xlabel('Quarter');
grid on













% graph for regression

% SVD regression


% Analyze the trend situation.
% figure
% detrend_dp=detrend(defaultFreq);
% trend = defaultFreq - detrend_dp;
% plot(defaultFreq);
% legend('Original Data','Location','northwest');
% hold on
% plot(trend,':r')
% plot(detrend_dp,'m')
% plot(zeros(size(defaultFreq)),':k')
% legend('Original Data','Trend','Detrended Data',...
%        'mean(Detrended)','Location','northwest')
% xlabel('Quarters'); 
% ylabel('DP');



% Test against the standard normal:chi2gof(SVD,'cdf',@normcdf)
% Use lillietest to determine if SVD follows a normal distribution: lillietest(SVD)
% The Lilliefors test is a 2-sided goodness-of-fit test suitable when a fully-specified null distribution is unknown and its parameters must be estimated. This is in contrast to the one-sample Kolmogorov-Smirnov test (see kstest), which requires that the null distribution be completely specified.





